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Abstract 

The recent fragmentation data for central collisions of Gold on Gold are even 
qualitatively different from those for central collisions of Argon on Scandium. 
The latter can be fitted with a lattice gas model calculation. Effort is made 
to understand why the model fails for Gold on Gold. The calculation suggests 
that the large Coulomb interaction which is operative for the larger system is 
responsible for this discrepancy. This is demonstrated by mapping the lattice 
gas model to a molecular dynamics calculation for disassembly. This mapping 
is quite faithful for Argon on Scandium but deviates strongly for Gold on 
Gold. The molecular dynamics calculation for disassembly reproduces the 
characteristics of the fragmentation data for both Gold on Gold and Argon 
on Scandium. 
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I. INTRODUCTION 



Recently we proposed a lattice gas model which was used to calculate mass distributions 
seen in heavy ion collisions at intermediate energy There are several features of this 

model that are attractive. The model can be used to study liquid-gas phase transition in 
mean-field theory and thus has links with Skyrme model studies of phase transitions. But 
the model can also be used to obtain cluster distributions and in this respect has close ties 
with percolation model of fragmentation which has been used with success in theoretical 
studies of heavy ion collisions. However, the lattice gas model has both kinetic energy and 
interactions, thus the scope of the model goes beyond that of standard percolation model. 

The model was used to fit the data obtained in Michigan State University for central 
collisions of Ar on Sc . There are several features of the model that are quite general and 
are also seen in experiments. First of all, there is a region of beam energy where the yield of 
Y{A) of fragments as a function of the mass number A of the fragment obeys a power law first 
noted in pioneering experiments by the Purdue group [0. This feature emerges in theoretical 
calculations also. In the lattice gas model the input in a calculation is the temperature and 
on general grounds the higher the beam energy, the higher the temperature. In the lattice 
gas model there is a temperature at which a power law will emerge. If in the vicinity of 
this critical temperature the mass distribution is fitted by a power law whose exponent is 
denoted by r then we expect to see a minimum in the extracted value of r at the critical 
temperature. It is difficult to miss a minimum in r in percolation type approach and also 
in the lattice gas approach. This minimum was seen in the Michigan experiment and was 
found in the calculation of at the experimental beam energy. Such a minimum in the 
value of r has also been seen in a recent experiment at Chalk River Nuclear Laboratories 
P] and the fit with a lattice gas model calculation is quite pleasing 0. Thus it seems that 
the lattice gas model simulates the fragmentation of nuclear matter reasonably well at least 
for medium mass collisions. 



Very recent data for Au on Au central collisions ||T0| are at variance with these general 
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expectations. In the range of beam energy 35 MeV/nucleon to 100 MeV/nucleon no min- 
imum in r was seen. Since this is now a much bigger system and phase transition effects 
should become more pronounced as we go to bigger systems, this absence of a minimum in 
r seems to cast serious doubts about the vahdity of the modeL The other remarkable result 
is that at 35 MeV/nucleon of beam energy the deduced value of r is much below 2; it is 
1.25. Such a low value of r will be difficult to obtain in the lattice gas model even including 
shape effects as discussed in ||rT|. Thus if we apply the lattice gas model as proposed in ^ 



and the calculation will not fit the data for the central collisions of Au on Au. 

The present work started as an effort to understand this puzzle. Between Ar on Sc 
and Au on Au the main difference is not only the sizes but also that in the latter case 
there is a huge Coulomb field whereas presumably in the former case the Coulomb field is 
merely a minor perturbation. Now the lattice gas model is a model of nearest neighbour 
interactions. The Coulomb field is a long range force and thus is not amenable to lattice 
gas type of approximation. One way to investigate the effect of the Coulomb field would 
be to try to add the effect of Coulomb force in the lattice gas type of approach but we 
could find no obvious way of doing that. We had provided a prescription, based on simple 
physical reasoning, to decide if two nucleons occupying neighbouring sites, form part of the 
same cluster or not. When in addition there is a long range force that distinguishes very 
strongly between neutrons and protons, this criterion clearly needs to be modified. There 
is no simple way of knowing how the modification should be made. Hence to find what 
effects a strong Coulomb field may have on a lattice gas model prediction, we need to try a 
somewhat convoluted approach. We first try to map the lattice gas model calculation to a 
molecular dynamics type calculation, both first done without any Coulomb interaction. If the 
calculations match quite faithfully then we can study the effects of the Coulomb interaction 
by adding that to the molecular dynamics calculation. Hopefully we will find that for Ar on 
Sc the Coulomb interaction does not severely change results and that the changes are very 
significant for the case Au on Au. We do not do an ab initio molecular dynamics calculation 
but only use it for disassembly from a thermally equilibrated source. The starting point 



of the lattice gas model is that matter has equilibrated at some temperature T. There 
are n nucleons and lattice sites (A^ > n). Which particular lattice sites are occupied 
are entirely dictated by statistical mechanics. Each cubic lattice has a size 1.0/po = 6.25 
fm^ and can, at most, be occupied by a single nucleon. This starting point of calculation 
will not be questioned in the present work but we will use two different prescriptions for 
obtaining the yield Y{A) against A. One of these two is our old prescription Nucleons 
occupying neighbouring sites will have attractive interactions — e and are considered to be 
part of the same cluster provided the kinetic energy of relative motion of these two nucleons 
does not overcome their binding. This is enough information to deduce the yield Y{A). In 
the alternative prescription that we will carry out here we fall back upon a more standard 
many body calculation. At the starting point when the nucleons have been initialised at 
their lattice sites and have their initial momenta, we will switch to a molecular dynamics 
calculation in which we will let the system evolve according to standard classical mechanics. 
Nucleons which stay together after arbitrarily long time are part of the same cluster. After 
a sufficiently long time the mass distribution is obtained. This can now be compared with 
lattice gas model predictions. Obviously for this test to be made with the lattice gas model 
we should choose for disassembly by molecular dynamics an interaction which suits the 
lattice gas model the best. The potential should be deepest with the value — e when the 
two nucleons are tq = 1.842 fm apart and should fall to before a/2 x 1.842 fm (the next 
nearest neighbour interaction is zero in the lattice gas model ). For distance less than 1.842 
fm the potential should quickly become repulsive ( two nucleons can not occupy the same 
site ). If the two prescriptions match for the yield Y{A), then we have linked the lattice gas 
model prediction for the yield Y{A) against A to a more well-known and better understood 
molecular dynamics approach. Considering that the lattice gas model can be easily linked 
with percolation model this in itself is quite interesting; we have provided a connection 
between percolation model results and molecular dynamics which seem to address totally 
different scenario to start out with. Secondly, in case the two results match in the absence of 
a Coulomb field, we can, in the molecular dynamics approach find out what a large Coulomb 
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field, which can not be incorporated in the lattice gas model, can do to the yield Y{A) since 
the Coulomb interaction is easily incorporated in molecular dynamics calculations. We will 
give the necessary details of these two calculations in the next section. 



We reiterate that our objective is not a molecular dynamics calculation as such 12,13 



we use disassembly by molecular dynamics with a nuclear force that produces results similar 
to those of lattice gas model. This work is focussed towards one question only: we ask why 
the lattice gas model which worked reasonably well for Ar on Sc fails in the case of Au on 
Au and if we can relate this failure to the large Coulomb field which is present in the second 
case. We note also that the Coulomb interaction can be incorporated in microcanonical 
models . 



II. DETAILS OF CALCULATION 

Motivation and details of the lattice gas model are given in [|I| and 0. For completeness 
some of these details are provided in this section. The calculations all require numerical 
simulation involving Monte-Carlo. The starting point of all our calculations is this. For a 
system of n nucleons we consider N lattice sites where N > n. is a parameter chosen in 
imj^ by requiring the best fit. The quantity N/n = Po/p where po is the normal density and 
p is the so called freeze-out density. It will be seen in the following that in the lattice gas 
calculation for fragments this freeze-out density is the actual density at which all clusters are 
calculated. For the molecular dynamics calculation that we will perform p is the density at 
which the initialisation is done according to prescription of equilibrium statistical mechanics. 
We then let the system evolve in time and the cluster distributions are calculated much later. 
Thus strictly speaking p is not a freeze-out density for molecular dynamics calculation but 
merely defines the starting point for time evolution. However since classical evolution of 
a many particle system is entirely deterministic, the initialisation does have in it all the 
information of the asymptotic cluster distribution. We will continue to call p the freeze-out 
density. 
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For initialisation, we assume that the nuclear part of the interaction is simply — e between 
nearest neighbours and zero otherwise. To begin a calculation we have to determine which 
of the sites nucleons occupy and what their momenta are. The two samplings can be done 
independently of each other. In a percolation model the N sites would be occupied with 
an occupation probablity p = n/N hj n nucleons in which each site has an equal apriori 
probability. Because of interactions this is somewhat more complicated in our case. Let us 
assume we are handling the case where we will take into account the Coulomb interaction 
explicitly in the initialisation. Starting with all lattice sites empty, the first nucleon ( a 
proton or a neutron as dictated by a Monte-Carlo decision ) is put at a site at random. If 
this first nucleon is a neutron then the probability of occupation of its nearest neighbours is 
proportional to exp(/?e) whereas all other sites have an occupation probability proportional 
to 1. As usual, P is the inverse of kT. These probabilities are now used to put in the second 
nucleon. If the first nucleon was a proton and the second one is a neutron then again the 
same probability of occupation will be used. But if the second nucleon is a proton also 
then the above occupation probabilities are changed to proportional to exp(/3e — (3uc) for 
the nearest neighbours and proportional to exp{—f3uc) for the other sites where Uc = e^/r, 
r being the appropriate distance between the two lattice sites. It is obvious how to repeat 
this procedure till the prescribed number of protons and neutrons are obtained. It is also 
obvious how to obtain the initial configuration when the Coulomb force is not explicitly 
included. In that case the prescription is identical with what was used in and [0. 

We do some initialisations where we take the Coulomb interaction explicitly and some 
initialisations when the Coulomb interaction is not taken into account separately. The 
nuclear part is always characterised by a strength — e which is the nearest neighbour type 
and has the same value irrespective of the isotopic spin. For a given nucleus the value of e is 
lower for the case when the Coulomb interaction is not explicitly added. This is required by 
demanding that the same binding energy is obtained in both the prescriptions. In the other 
case when the Coulomb interaction is explicitly included the nuclear part of the interaction 
will lead to a larger binding which will be somewhat compensated by the repulsive Coulomb 
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part. For the cases dealt with here e changes by roughly 1 MeV. 

For disassembly by molecular dynamics we approximate the nuclear part of the force by 
a well-known parametrisation ||16|| : 



Here tq is the distance between the centres of two adjacent lattices. We have chosen p = 2, 
q = 1 and a = 1.3. The other constants A and B are chosen so that the potential acquires 
the prescribed value — e at r = tq. With this potential the interaction between two nucleons 
is zero when they are more than 1.3ro apart and the interaction begins to become strongly 
repulsive for r significantly less than tq ; yet the potential is smooth enough that accurate 
numerical solutions of time evolution of nucleons can be obtained. The time evolution 
equations for each nucleon are, as usual, given by dpi/dt = —Y.jjf^i'^i'^ij'ij) and dri/dt = 



The lattice gas predictions for cluster production can only be calculated for the case 
where the Coulomb interaction is not explicitly included but only through a lower value of e. 
As mentioned already the cluster distribution is calculated immediately after initialisation. 
The lattice filling is done and the momenta are then generated from a Monte-Carlo sampling 
of a Maxwell-Boltzmann distribution with a prescribed temperature. Nucleons in two neigh- 
bouring cells are considered to be part of the same cluster if the kinetic energy of relative 
motion is not large enough to overcome the attractive interaction, i. e., — e < 0. 

Here fi is the reduced mass and is equal to m/2. This definition is the simplest that one 
can provide and is physically reasonable. It should be emphasized that it is by no means a 
unique one. The prescription manages to reduce a many body problem of cluster production 
into a sum of independent two body problems. One can easily construct scenarios where 
this prescription may underestimate the size of a cluster and scenarios where this prescrip- 
tion may overestimate the size of a cluster. It should be pointed out that this formula for 
bond formation has the same structure as the one used in [Q. Since each particle obeys 
the Maxwell-Boltzmann distribution, the distribution of relative momentum between two 
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A[B{rQ/ryP — (ro/r)''] exp[l/(r/ro — a)], for r /tq < a 



for r/ro > a 



(1) 



Vi/m. 



particles is also a Maxwell-Boltzmann, i.e., P(pr-) = [l/(27r/iA;T)^/^] exp[— We 
can then write down a formula for the bonding probability which is temperature dependent 

Switching to a variable E = pl/2n we get 

which is identical with the formula of 0. 

For molecular dynamics calculation after initialisation we do the time propagation long 
enough so that for cluster production the asymptotic time has been reached. Two nucleons 
are part of the same cluster if the configuration distance between them is less than 1.3ro. 
We stop the calculation after the original blob of matter has expanded to 64 times it volume 
at initialisation. For low temperature this means doing the time evolution as long as 1000 
fm/c. We use a time step of 0.1 fm/c and update positions and momenta half a time step 
apart ("leap frog" method). The energy conservation in our calculation is accurate to within 
1 one percent. The program conserves total momentum identically. Below we now consider 
specific cases. 



III. RESULTS 

In Ref. [0] we found that a freeze-out density p = .39po gave the best fit with data. 
Here we present data with this freeze-out density. A few calculations with a higher value of 
freeze-out density were also performed but only to ascertain that the trends of the results 
are not strongly dependent on the freeze-out density employed. 

Fig. 1 shows the results of a Y{A)-A plot of a lattice gas calculation. This is a repeat 
of the type of calculation done in |jl]J^. The value of e used is 3.7 MeV. This curve should 
be compared with a molecular dynamics calculation shown in Fig. 2. This calculation uses 
the same e and no explicit Coulomb interaction. The similarity between Figs. 1 and 2 is 
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quite striking and leads us to conclude that the simple prescription of cluster counting is very 
reasonable. In Fig. 3 we have done a molecular dynamics calculation where we explicitly put 
in the Coulomb interaction. Accordingly the value of e has been increased from 3.7 MeV to 
4.7 MeV. There are now some changes from the results of Figs. 1 and 2, but not a great deal. 
Specially the deduced value of the slope r is again the lowest at T = Tc and rises both below 
and above this temperature. Tc is 1.1275e and is the critical temperature in the lattice gas 
model. The results for r are summarised in Fig. 4 where it is seen that the explicit inclusion 
of the Coulomb interaction has not modified the predominant characteristics observed in 
calculations without explicit inclusion of the Coulomb interaction. Here r was obtained 
from hnear fits of fragmentation distributions in \ogY{A) vs log A plots. That is, r is 
determined by minimizing the defined as: 

X' = E(^(A) - F,f. (4) 

i 

Here F(Ai) = \ogY(Ai) — const — rlogAi is the fitted yield for fragment of size Ai, and 
Fj = logFj is the corresponding simulated yield. To maintain sufficient statistics and to 
exclude the largest cluster, only fragments of sizes between 1 and 12 were used. In Fig. 5 
for completeness we have shown aY{Z) against Z curve. This is the type of curve that are 
typically presented as experimental results. 

Figs 6, 7 and 8 give our results for Au on Au. We plot both Y{A) against A and Y{Z) 
against Z. For molecular dynamics without explicit inclusion of the Coulomb interaction we 
have used e = 3.7 MeV and for calculation with explicit inclusion of the Coulomb interaction 
we have used e = 4.7 MeV. However now the results are very different for the two cases. In 
one (Fig. 6) there is a minimum in r at T = Tc. A second spike at T = OATc is indicative 
of a percolating cluster. In Fig. 7 with explicit inclusion of the Coulomb interaction, the 
percolating cluster has disappeared. We also see that below the critical temperature the r 
values from the two calculations begin to diverge. With Coulomb explicitly included the 
minimum in r has disappeared and one can get a value of r much less than 2. These results 
are summarised in Fig. 8. Our results are in qualitative agreement with the experimental 
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results for Au on Au. Experimental results are given as a function of the beam energy and 
thus we need a conversion from temperature to beam energy. This conversion needs to be 
done carefully because at initialisation, which is the starting point of our calculation we 
can expect that some energy is already in collective motion and does not appear as thermal 
excitation. The model does not include this aspect. In Ref. a phenomenlogical mapping 
of temperature to beam energy was deduced from [^. As an estimate only if we assume 
that at the initial time 3/8 of the initial energy is stored in collective motion then a beam 
energy of 35MeV/nucleon would correspond to O.STc. Our calculated value of r is then 1.4 
compared to the experimental value of 1.25 |jlO|- As in experimental data the calculated 
Y{A) against A deviates from a power law with higher excitation energy. We nonetheless 
deduce a effective value of r from a very approximate fit and these are shown in Fig. 8. We 
regard r as a measure of global feature of Y{A) , although power-law fits are poor at high 
temperatures. To maintain sufficient statistics, we used fragments of size between 1 and 
20 for high temperatures (T > Tc) when heavier fragments are rare. For low temperatures 
(T < Tc) larger fragments were also included. Our calculation at about I.IT^ fits the data for 
beam energy 100 MeV/nucleon. For the calculation with Coulomb interaction included we 
use Tc merely as an energy scale; there is no implication that Tc is the critical temperature 
of the system. The principal point we want to emphasize is that we have reproduced the 
most significant features of the data for Au on Au as contrasted with those for Ar on Sc, 
namely that in the former case there is no minimum in the value of r and that the value of 
r can be significantly below 2. 



IV. CONCLUSION 

This problem started out as an effort to understand why the fragmentation data for Au 
on Au are so different from that of Ar on Sc and if the data totally ruin all validity of the 
simple concepts used in the lattice gas model. The calculation done here suggests that the 
lattice gas model is reasonable for medium mass collisions; it probably would have been as 
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valid for collisions of very large masses but for the very large Coulomb force which begins to 
make its presence felt and destroys the simple predictions. It has often been assumed that 
the larger the system of colliding nuclei the better is the chance of learning about phase 
transition in nuclear matter. However larger colliding masses also bring in much larger 
Coulomb forces and it will be necessary to take into account of the Coulomb effects before 
the signals for phase transitions can be understood. With large masses the mean field of 
the protons are very different from that of the neutrons and theories must be able to treat 
them differentially. There clearly are needs for simple theories which are able to handle this 
difference. 
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FIGURE CAPTIONS 



Fig. 1 The mass yield distributions obtained from the lattice gas model for lattice N — 6^ 
and n = 85, at temperatures T/Tc =0.5, 1.0 and 1.5. Here Tc = 1.1275e is the thermal 
critical temperature. 

Fig. 2 The mass yield distributions obtained from molecular dynamics calculations with- 
out the inclusion of Coulomb interaction. The lattice size, number of nucleons and 
temperatures are the same for the corresponding curves in Fig. 1. 

Fig. 3 The same as Fig. 2, but the Coulomb interaction is taken into account. 

Fig. 4 The value of r obtained from lattice gas model and molecular dynamics calculations 
are plotted as a function of temperature for N — 6^ and n = 85. 

Fig. 5 The charge yield distributions obtained from the molecular dynamics calculation 
with Coulomb interaction are shown. 

Fig. 6 The mass and charge yield distributions for Au on Au collisions obtained from the 
molecular dynamics without Coulomb interaction. 

Fig. 7 The same as Fig. 6, but with Coulomb interaction. 

Fig. 8 The value of r in Au on Au collisions obtained from molecular dynamics calculations 
with and without Coulomb interaction are plotted as a function of temperature. 
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